Statistical features of freely decaying two-dimensional hydrodynamic turbulence 
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Statistical characteristics of freely decaying two-dimensional hydrodynamic turbulence at high 
Reynolds numbers are numerically studied. In particular, numerical experiments (with resolution 
up to 8192 x 8192) provide a Kraichnan-type turbulence spectrum Ek ~ k~ 3 . By means of spatial 
filtration, it is found that the main contribution to the spectrum comes from the sharp vorticity 
gradients in the form of quasi-shocks. Such quasi-singularities are responsible for a strong angular 
dependence of the spectrum owing to well- localized (in terms of the angle) jets with minor and/or 
large overlapping. In each jet, the spectrum decreases as A: -3 . The behavior of the third-order 
structure function accurately agrees with Kraichnan direct cascade concept corresponding to a 
constant enstrophy flux. It is shown that the power law exponents ( n for higher structure functions 
grow more slowly than the linear dependence of n, which testifies to turbulence inter mittency. 
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1. Introduction. It is known that developed two- 
dimensional hydrodynamic turbulence, in contrast to 
three-dimensional turbulence, in the inert ial interval 
of scales possesses an additional integral of motion — 
enstrophy. Enstrophy is an integral of the squared 
vorticity f Q 2 dr. As Kraichnan demonstrated in 
1967 [1 j, the existence of this integral generates in 
the inertial interval its own Kolmogorov spectrum of 
turbulence: 

E(k) ~ k~ 3 , (1) 

(now called the Kraichnan spectrum). This spec- 
trum corresponds to a constant enstrophy flux to- 
ward the region of small scales. Simultaneously, ac- 
cording to pp, a conventional Kolmogorov spectrum 
E(k) ~ /c -5 / 3 with a constant energy flux directed 
toward the region of small values of k (inverse cas- 
cade) is formed in the region of large scales. After 
that paper was published, many numerical experi- 
ments were performed (we note only some of them: 
[2]- [9]; a more complete list can be found in [TU]), 
which testified in favor of existence of the Kraichnan 
spectrum. At the same time, already in the first nu- 
merical experiments there was observed (see, e.g., [2]) 
the emergence of sharp vorticity gradients (consistent 
with the high Reynolds number). It corresponds to 
the formation of jumps similar to shock waves which 
thickness is small as compared with their length. 
Based on this observation, Saffman [11 proposed an- 
other spectrum, which differs from the Kraichnan dis- 
tribution: E(k) ~ k~ A . The Saffman spectrum was 
obtained under the assumption that the main con- 
tribution to the spectrum comes from isotropically 



distributed vorticity shocks. On the other hand, it 
follows from simple considerations that the spectrum 
with such shocks is expected to have the Kraichnan 
type behavior rather than the Saffman distribution. 
Indeed, if one considers the Fourier transform of a 
vorticity jump ft oc 0(x), where 0{x) is the Heaviside 
function, then one can easily obtain that 
immediately yields the spectrum E(k) ~ k~ 3 . How- 
ever, the situation is not too simple. If one assumes 
that the characteristic length of the step L ^> 
then the energy distribution from one such step in 
the /c-space has the form of a jet with an apex angle 
of the order of (ZcL) -1 . The maximum value of the 
energy distribution along the jet decreases as ex k~ 3 
[lOj [I2j [13] . As was demonstrated in [lOj [12] , such an 
estimate at kL ^> 1 follows from the rigorous analysis. 

The main attention in this work is paid to the nu- 
merical study of statistical characteristics of the di- 
rect cascade of freely decaying two-dimensional tur- 
bulence and the influence of vorticity quasi-shocks on 
these characteristics. As compared to previous works, 

(i) the spatial resolution is significantly increased 
(up to K x K = 8192 2 points). High-resolution 
numerical simulations confirmed all previous results 
[IS [13]; 

(ii) it is proved numerically that the spectrum ([!]) 
arises owing to vorticity quasi-shocks, which form a 
system of jets with weak and strong overlapping in 
the /c-space; 

(iii) in each of these jets, the decrease in E(k) at 
high values of k is proportional to A: -3 , which yields 
the spectrum ([I]) after averaging over the angle; 

(iv) the third-order structure function ( (<5^) 2 <5^|| ) 
depends linearly on R = |r — r'| (where 5Q = Q(r) — 
^(r'), 5v\\ is the projection of the velocity difference 
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onto the vector R), in complete agreement with the 
Kraichnan theory. The higher structure functions 
of velocity have the exponents £ n , which grow more 
slowly than the first power of n, testifying to turbu- 
lence intermittency. 

2. Main equations and numerical scheme. 

As was noted in Introduction, the emergence of sharp 
gradients of vorticity was observed in the majority 
of numerical experiments (see, e.g., [2]- [9]) modeling 
two-dimensional turbulence at high Reynolds num- 
bers, i.e., in the regime where the Euler equation can 
be considered in the zero approximation instead of 
the Navier-Stokes equation. 

In terms of vorticity £1 = dv y /dx — dv x /dy, the 
Euler equation is written as 



on 

~dt 



(vV)ft with div v = 0, 



(2) 



i.e., Q is a Lagrangian invariant. As was demon- 
strated in [TO] [12], the reason for the emergence of 
strong gradients of vorticity in two-dimensional tur- 
bulence is connected with compressibility of the field 
B = rot (fii), which was called divorticity. This vec- 
tor is tangent to the level lines fi(r) = const, and 
the field B is frozen into the fluid, i.e., it satisfies the 
equation [6] [14] 



dB 



= rot [v x B] . 



(3) 



It follows from Q that B changes owing to the veloc- 
ity component v n normal to B, which is the velocity 
of the lines of the field B due to its frozenness. The 
Lagrangian trajectories determined by v n , 



dr 

dt 



v n (r,t); r| £ = = a, 



define the mapping 

r r(a,t). 

In this case Eq. (|3| admits integration [10] [12] : 
(B (a)- Va)r(a,f) 



B(r,f) 



J 



(4) 



(5) 



(6) 



where B (a) is the initial value of the field B, and J 
is the Jacobian of mapping (J5|. As the trajectories 
([5| are defined by the normal component of velocity 
rather than by velocity proper, the Jacobian J is not 
fixed. It can take arbitrary values, in particular, it 
can tend to zero. This implies compressibility of the 
field B. It should be noted that Eq. (|6| is an analog of 
the vortex lines representation, which was introduced 
in p~5j [16] for the three-dimensional Euler equation. 



It is known that the emergence of discontinuities 
in gas dynamics is associated with vanishing of the 
Jacobian of the corresponding compressible mapping, 
i.e., the transformation from the Eulerian to the La- 
grangian description. Just the decrease in J is respon- 
sible for the occurrence of sharp gradients of vorticity 
resulting from compressibility v n : divv n ^ 0. In the 
two-dimensional incompressible Euler hydrodynam- 
ics, however, singularity formation is impossible in 
a finite time [17] ; possibly, it can occur in an infinite 
time. Numerical experiments [10] [13] yield an expo- 
nential increase in vorticity gradients with the charac- 
teristic increment of the order of the maximum value 
of vorticity rather than double-exponential growth, 
as it follows from the estimates made in [18 . The 
same behavior is observed in our numerical experi- 
ments, which ensure a better resolution (see the next 
paragraph) . 

Let us briefly describe the numerical scheme used. 
The equation for vorticity ([2| was solved numerically 
in a square box with periodic boundary conditions 
along both coordinates. The box size was chosen to be 
equal to unity. To eliminate bottleneck instability and 
ensure a sink of energy at small scales, we introduced 
dissipation into the right hand side of Eq. ([2| , which 
was taken in the form of hyperviscosity 

(_!)«+ V„V 2 "c, » n = 10- 20 (^) 2 " , n = 3. 

The use of hyperviscosity allows appreciable reduc- 
tion of the range of scales where the influence of dis- 
sipation is significant and, therefore, expansion of the 
inertial region of the spectrum [T5] . When the spatial 
resolution was changed, the hyperviscosity coefficient 
was scaled in such a manner that the value of dissi- 
pation for the shortest wave perturbations that could 
be presented on the computational grid remained con- 
stant. The effect of hyperviscosity on the flow as a 
whole always remained extremely small: in the nu- 
merical experiments, the total energy was preserved 
within 0.01 % . 

Equation ^ with hyperviscosity was solved by the 
pseudo-spectral Fourier method [20], i.e., the deriva- 
tives were calculated and the Poisson equation (for 
finding the stream function from vorticity) was solved 
in the spectral space, whereas the nonlinear terms 
were calculated on a computational grid in the phys- 
ical space. The transformations from the physical to 
the spectral space and back were performed through 
the Fast Fourier Transform (FFT) with the FFTW 
library used for that purpose [2T] , The use of hy- 
perviscosity gives rise to a stiff term in the equation; 
because of this term, the time step of integration is 
restricted for any explicit scheme by a prohibitively 



FIG. 1: Initial distribution of vorticity, N = 20 



FIG. 2: Distribution of vorticity at t = 100, N = 20 



small value (~ l/if 2n ). To avoid it, time integration 
was performed by a hybrid third-order Runge-Kutta 
/ Crank-Nicholson method [22] . The convective term 
is approximated explicitly, and an implicit scheme is 
used for the dissipative (linear) term. 

The computations were performed both on a con- 
ventional multiprocessor cluster (for which the code 
was parallelized with the help of the MPI library; up 
to 128 processors were used) and on a GPU cluster 
(with the use of the NVIDIA CUDA technology) at 
the Computer Center of the Novosibirsk State Uni- 
versity. The spatial resolution was up to 8192x8192 
points. 

3. Numerical experiments. The initial condi- 
tion was chosen in the form of two sets of vortices 
with positive and negative vorticity. All vortices had 
a Gaussian shape with the maximum value of |0| 
equal to unity and with random sizes (the vortex radii 
were uniformly distributed in a certain interval); the 
vortex centers were also randomly distributed over 
the entire domain. The total vorticity was equal to 
zero. The number of positive vortices was equal to 
the number of negative vortices; this number varied 
depending on the spatial resolution from N = 10 at 
2048x2048 to N = 40 at the maximum spatial res- 
olution. Variations of the initial conditions did not 
change the qualitative behavior of the vortices and 
their statistical characteristics. A typical initial con- 
dition is shown in Fig. 1. 

Figure 2 shows the vorticity distribution at t = 



100. This time is close to the maximum value t max 
at which the computation was terminated. The time 
t = tmax is determined from the condition that the 
edge of the energy spectrum kb oun d reaches 2/3 of the 
maximum value k max in the box. At t > t max , the 
aliasing phenomenon arises [20], where the spectrum 
is distorted because harmonics located at a distance 
of 2k max from each other in the spectral space cannot 
be distinguished on the computational grid. 

As is seen from Fig. 2 (resolution 8192 2 ), the struc- 
ture of vortices at t = 100 is substantially changed ow- 
ing to interaction of vortices, emergence of shear flows 
generated by vortices, and their stretching. Figure 3 
shows the angularly averaged compensated spectrum 
of turbulence k 3 E(k) at different instants of time. At 
times close to t max: a plateau is formed in the com- 
pensated spectrum, i.e., a Kraichnan-type spectrum 
E(k) ~ k~ 3 is formed almost on two decades. 

As was noted previously (see also jTQl [13]), such a 
spectrum is formed owing to the sharp vorticity gra- 
dients. Figure 4 shows the distribution of |B| = 
at t = 100, which demonstrates a significant increase 
(~ 100) in the initial gradients in the vicinity of lines 
corresponding to the positions of vorticity shocks. Be- 
tween the lines, |B| practically does not increase. At 
greater times (t > 100), the lines become consid- 
erably longer; correspondingly, the pattern contains 
more lines, and the distribution of these lines becomes 
rather twisted (turbulent). 

At the initial stage, at several turnover times n 
(T* = n27r/Q max ), |B| grows almost exponentially, 
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and then the saturation is observed. Thus, for Fig. 
5, the number of turnovers is n ~ 4, and |B| increases 
almost by a factor of 300. 

A significant increase in vorticity gradients and in 
the density of lines corresponding to the positions 
of vorticity shocks testifies that their role in spec- 
trum formation is rather important. This (princi- 
pal) issue was not completely clarified in [lOj H3] , 
To verify numerically that the Kraichnan-type spec- 
trum E(k) ~ k~ 3 is generated by sharp gradients 
of vorticity, we performed spatial filtration, in con- 
trast to [TUJ [13] , where filtration in the /c-space was 
performed. In the distribution of the field B(r), we 
retained only those regions where |B| exceeded some 
threshold value Bo > 0. In regions with |B| < £?o, 
the field B was set to zero. Using the "filtered" field 
B (r), we found the Fourier transform B&. Further 
the Fourier transform of the filtered velocity was de- 
termined: 



1 



kx 



kxB* 



and then the "filtered" spectrum E(k). (It is evident 
that this procedure introduces its own contribution 
to the spectrum E(k); this contribution is ensured 
by the shocks |B| generated by power-law tails de- 
creasing at least as k~ A . In other words, filtration 
alters only the long- wave part of the spectrum.) 
Figure 6 shows the "filtered" spectra k 3 E{k) for dif- 
ferent threshold values Bo. As a result_of such fil- 
tration, the (compensated) spectrum k 3 E(k) experi- 
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FIG. 5: Maximum value of |B| versus time for the initial 
condition shown in Fig. 1 (logarithmic scale, the straight 
line corresponds to the exponential growth) 



enced changes only in the region of small values of k 
and remained unchanged in the region of large values 
of /c, which is a numerical consequence of the fact that 
the Kraichnan-type spectra are induced by vorticity 
quasi-singularities. 

The reason for the emergence of the dependence 
k~ 3 for the spectrum is closely connected with its an- 
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FIG. 6: Filtered compensated spectra k 3 E(k) for different 
threshold values i?o 



FIG. 7: Distribution of the two-dimensional spectrum 



gular distribution. Figure 7 shows the distribution of 
the two-dimensional spectrum e(k x ,k y ). This distri- 
bution considerably varies with the angle. It is a num- 
ber of jets with weak and/or strong overlapping. This 
behavior is also clearly visible in the angular depen- 
dence of the spectra with different values of k (Fig. 8). 
Along each jet, the spectrum e& changes as /c -4 , which 
yields (after angular averaging) a Kraichnan-type dis- 
tribution E[k) 0. 

Apart from the spectrum, important characteristics 
of turbulence are structure functions. As was noted 
above, at times t close to t max , the distribution of the 
lines of the maximum values of |B| becomes substan- 
tially more complicated; we can say that it becomes 
turbulized. In particular, the angular overlapping of 
jets becomes more noticeable, and the spectrum be- 
comes more isotropic. To elucidate whether the tur- 
bulence considered here is close to the Kraichnan-type 
turbulence, we performed numerical experiments at t 
close to t m ax an d obtained the dependence of the cor- 
relation function 

D(R) = ((fi(r) - fi(r')) 2 Qv(r) - v(r')] • [R]) /R) 

on R = |r — r'|. For the direct cascade, this depen- 
dence on i?, as was predicted by Kraichnan (see, e.g., 
review [23]), should be linear with a proportionality 
coefficient 77, which is the flux of enstrophy toward 
the region of small scales. Figure 9 shows that this 
dependence on R is indeed close to linear, especially 
at small values of R. At large values of R, there is a 
small inflection, which, in our opinion, is associated 



with deviation of the spectrum from the isotropic one, 
i.e., with the jet-like structure of the spectrum. 

In addition to the correlation function D(R), we 
also calculated the structural functions of velocity 

5„(i?) = ([(v(r')-v(r))-^]"). 

It is known that the dependence of the power of the 
structure functions ( n (S n (R) ~ R^ n ) on n character- 
izes the degree of turbulence intermittency. For the 
Kolmogorov-Kraichnan regime, ( n is a linear depen- 
dence of n. Our numerical experiment shows that the 
local dependence ( n (Fig. 10) deviates from a linear 
dependence and has a tendency to saturation, espe- 
cially at large values of i?, which testifies to intermit- 
tency of this type of turbulence. It should be noted 
that all ("n m the absence of the dissipative scale de- 
pend linearly on n as R — ^ (see Fig. 10). It is caused 
by the analytical dependence of velocity difference at 
small scales. 

4. Conclusion. The main resume of this work is 
the fact that the Kraichnan-type power-law spectrum 
is formed due to quasi-singularities, which appear in 
solving the Cauchy problem for the two-dimensional 
Euler equation. Though the collapse as the pro- 
cess of singularity formation in a finite time is for- 
bidden in two-dimensional hydrodynamics, but there 
is a tendency to the emergence of such singularities. 
As it follows from our numerical experiments, vor- 
ticity gradients increase by more than two orders of 
magnitude. In our opinion, this process can be en- 
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FIG. 8: Compensated spectrum k A e versus the angle (j) for ko = 500, ko = 1000, and ko = 1500 
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FIG. 9: Correlation function D(R) 



FIG. 10: Exponents £ n (local) as functions of R 



hanced in the presence of forcing. This assumption 
is based on the presence of a logarithmic correction 
in terms of k in the Kraichnan spectrum, which is 
associated with a weak interaction nonlocality. It is 
also of interest that the third-order correlation func- 
tion D(R) demonstrates the Kraichnan- type behavior 
at the stage of turbulence decay, when a plateau is 
formed in the compensated spectrum k 3 E(k), which 
does not contradict the conclusion that spectrum ([I]) 
is generated by vorticity quasi-shocks. At a suffi- 
ciently late stage of turbulence development, however, 
the turbulence spectrum has a rather anisotropic dis- 
tribution over the angle owing to jets. 
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